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SECTION  I 


INTRODUCTION 

In  studying  foreign  object  damage  (FOD)  of  laminated  fiber-reinforced 
composites,  one  often  encounters  extremely  complex  material  as  well  as 
structural  responses.  First  of  all,  the  high  velocity  impact  usually 
results  in  a substantial  permanent  deformation  in  the  contact  region. 
Secondly,  the  structural  response  is  transient  and  it  often  requires  a 
large  number  of  vibrational  modes  to  ensure  an  adequate  description  of  the 
motion  [1].  Consequently,  most  of  the  analytical  work  [2-3]  is  restricted 
to  the  elastic  response  of  the  laminate  by  neglecting  the  plastic  deformation 
and  by  assuming  that  the  contact  force  between  the  projectile  and  the 
laminate  is  known.  Even  with  such  simplifying  assumptions,  none  is  able  to 
predict  quantitatively  the  damage  of  the  laminate  due  to  impact. 

In  contrast  to  the  small  quantity  of  analytical  work,  experimental 
work  is  abundant  [4].  Various  tests  have  been  performed  in  an  effort  to 
understand  the  failure  mechanisms  with  the  objective  of  developing  practical 
protective  measures.  Although  testing  yields  quantitative  results,  it 
is  In  general  very  restrictive  and  expensive.  In  view  of  the  great 
variety  of  composites,  a reliable  analytical  method  is  highly  desirable 
in  general  applications,  especially  In  optimal  design. 

The  present  work  alms  at  developing  a finite  element  method  to 
determine  the  dynamic  response  as  well  as  the  energy  that  causes  damages 
in  the  laminated  composite  due  to  Impact  of  foreign  objects.  Such  damage 
energy  can  be  used  to  estimate  the  residual  strength  of  the  composite 
after  Impact.  In  this  report,  a higher  order  beam  finite  element  is 
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presented  complete  with  stiffness  and  consistent  mass  matrices. 

The  finite  element  is  evaluated  by  considering  homogeneous  and  iso- 
tropic beams  subjected  to  impact  of  elastic  spheres  for  which  ana- 
lytical solutions  are  available.  The  elements  are  then  employed 
to  study  beams  of  a glass-epoxy  laminated  composite.  Solutions 
are  compared  with  the  experimental  data  obtained  by  Husman  et  al  [5]. 


2 


■ r 


SECTION  II 


THE  FINITE  ELEMENT  METHOD 

We  will  consider  symmetrically  laminated  cross-ply  composites. 
Extension  to  unsymmetric  laminates  which  exhibit  bending-extension  coupling 
can  be  easily  carried  out  in  the  same  manner.  In  the  following,  we  will 
derive  the  stiffness  matrix  and  the  mass  matrix  for  a homogeneous  and 
isotropic  beam  element  for  which  evaluation  of  the  element  efficiency  will 
be  made.  The  finite  element  for  a symmetrically  laminated  cross-ply 
composite  can  be  obtained  from  the  isotropic  one  by  using  the  appropriate 
stiffness  coefficients.  Bernoulli -Euler  beam  theory  will  be  employed 
in  modeling  the  finite  element. 

1.  Finite  Element  Model  for  Isotropic  Materials 

The  following  power  series  for  the  shape  function  for  the  beam 
element  in  flexural  deformation  is  taken: 

v = a-j  + a2x  + a^x^  + a^x^  + a5x4  + a^x5  (1) 

With  this  displacement  function,  there  are  three  degrees  of  freedom  at 
each  nodal  point,  namely,  the  transverse  displacement  v^ , the  rotation 
ei  and  the  curvature  (1=1,2).  The  displacement  function  can  be 
expressed  in  terms  of  the  six  generalized  nodal  displacements. 

Consider  an  elastic  beam  element  of  length  L,  cross-sectional  area 
A,  moment  of  inertia  I and  the  elastic  modulus  given  by  E^.  The  potential 
energy  U and  the  kinetic  energy  T stored  in  this  element  are  given  by 


(2) 


and 


L 

(tf)2dx 

0 


(3) 


respectively.  In  Eq.  (3),  p is  the  mass  density,  and  a dot  represents 
the  derivative  with  respect  to  time. 

Corresponding  to  each  degree  of  freedom  at  the  nodal  point,  there  is 
a generalized  force.  In  this  case  we  have  , m.  and  corresponding  to 
v.,  e..  and  K.j , respectively.  The  quantities  Q.  and  m.  are  the  usual  shear 
force  and  moment,  and  p^  is  a hypermoment  with  the  dimension  lb-ft2. 

The  potential  energy  of  the  generalized  end  forces  for  the  element 
is 

2 

We  = " £ (Vi  + miei  + piV  (4) 


Noting  that  the  displacement  function  of  the  element  is  now 
expressed  in  terms  of  the  six  nodal  displacement  components,  we  can  write 
U and  T given  by  Eqs.  (2)  and  (3)  in  terms  of  the  six  nodal  displacement 
components.  The  equation  of  motion  can  be  obtained  by  applying  Hamilton's 
principle,  i.e.. 
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(U-T+Wg)dt=0 


(5) 
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with  <$a.  = Oat  tQ  and  t^ . 

The  result  is  a system  of  six  equations  which  can  be  written  in 
matrix  form  as 


{F}  = [k]  {A}  + [m]  {A} 


(6) 


where  { F > is  the  vector  of  the  generalized  forces  associated  with  the 
nodal  degrees  of  freedom  {Al;[k]  is  the  element  stiffness  matrix  and  [m] 
is  the  element  mass  matrix.  The  explicit  expression  of  Eq.  (6)  is  given  by 
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From  the  d'Alembert's  dynamic  condition  of  equilibrium  and  the 
condition  of  continuity  at  the  nodal  points,  we  obtain  the  system  of 
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MW 


- - - a mmam 


(12) 


Muw  = (PW  - t«lblt 


where 

[H]  = [MJ  + [K]  (13) 

and 

{blt  = {&)t  + At{A)t  + {&}t  (14) 

The  solution  for  l4lt+At  can  be  obtained  by  inverting  [H]  in  Eq.  (12). 
However,  a more  efficient  scheme  is  to  solve  Eq.  (12)  by  the  Gauss 
elimination  procedure.  With  Eqs.  (9)  and  (10),  the  current  generalized 
displacement,  velocity,  and  acceleration  are  obtained  from  those  at  the 
previous  time.  If  appropriate  initial  conditions  are  given,  the  total 
response  history  can  be  generated  by  this  step  by  step  procedure. 

We  have  also  employed  another  finite  difference  procedure  suggested 
by  de  Vogelaere  [8]  in  some  of  the  evaluative  examples  to  be  discussed 
later,  and  compared  with  the  present  method.  It  was  found  that  these  two 
methods  were  comparable. 

It  should  be  noted  here  that  the  choice  of  the  time  Increment  At  Is 
very  critical  In  assuring  the  stability  of  the  solution.  According  to 
the  conclusion  of  Leech,  Hsu  and  Mach  [9],  solution  for  Eq.  (8)  would 
be  stable  If 

At  2 2/u>m 


(15) 


equations  of  motion  for  the  whole  structure: 


(P)  = [K]{A)  + [M]{A} 


(8) 


where  { P } indicates  the  external  loads,  and  [K]  and  [M]  are  the  assembled 
structural  stiffness  and  mass  matrices,  respectively. 

2.  Integration  of  Time  Variable 

The  solution  for  the  equations  of  motion  given  by  Eq.  (8)  can  be 

solved  by  various  methods  [6].  A simple  finite  difference  form  proposed 

by  Wilson  and  Clough  [7]  will  be  followed  in  the  present  analysis. 

The  essence  of  the  procedure  proposed  by  Ref.  [7]  resides  in  the 

assumption  that  the  variation  of  the  acceleration  (A)  Is  linear  in  the 
time  interval  At.  Using  the  assumption,  we  have 


{A} 


t+At 


{Alt  + $ 


(A}4 


At 

T 


(A) 


t+At 


(9) 


and 


{A}t+At  = {A}t  + At{A}t  + {A}t  + T-  {A}t+At  (10) 


Substitution  of  equations  (9)  and  (10)  in  (8)  yields 
■ M«»)t  * 4t{i)t  * <4,t  * 


(11) 


from  which  we  obtain 
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= {PW  " ^Kj{b}t 

where 

[H]  = [M]  + tK]  (13) 

and 

{ b > t a { A} t + At{A>t  + { A)t  (14) 

The  solution  for  can  be  obtained  by  inverting  [H]  in  Eq.  (12). 

However,  a more  efficient  scheme  is  to  solve  Eq.  (12)  by  the  Gauss 
elimination  procedure.  With  Eqs.  (9)  and  (10),  the  current  generalized 
displacement,  velocity,  and  acceleration  are  obtained  from  those  at  the 
previous  time.  If  appropriate  initial  conditions  are  given,  the  total 
response  history  can  be  generated  by  this  step  :y  step  procedure. 

We  have  also  employed  another  finite  difference  procedure  suggested 
by  de  Vogelaere  [8]  in  some  of  the  evaluative  examples  to  be  discussed 
later,  and  compared  with  the  present  method.  It  was  found  that  these  two 
methods  were  comparable. 

It  should  be  noted  here  that  the  choice  of  the  time  Increment  At  is 
very  critical  in  assuring  the  stability  of  the  solution.  According  to 
the  conclusion  of  Leech,  Hsu  and  Mach  [9],  solution  for  Eq.  (8)  would 
be  stable  if 

At  ^ 2/%  (15) 

where  is  the  highest  natural  frequency  of  the  beam  element.  They 
also  suggested  that  At  should  be  no  greater  than  one-sixth  of  the  critical 

7 


r 

i 


SH1’ 


value  in  order  to  obtain  the  fine  high-frequency  detail  of  the  response. 
However,  we  have  found  that  to  insure  the  stability  and  to  contain  the 
higher  frequencies  of  the  higner  order  beam  element  model  At  should  be 
no  greater  than  one-tenth  of  the  critical  value.  For  nonlinear  problem 
At  should  be  even  smaller. 

3.  Evaluative  Example  - Impulsive  Loading 

In  order  to  test  the  efficiency  of  the  present  higher  order  element 
we  consider  a simply-supported  beam  of  length  L subjected  to  a sine  pulse 
given  by 

P * P0  sin  for  0 - t - t (16) 

P * 0 for  t - t (17) 


The  force  is  located  at  the  center  of  the  span.  An  analytical  solution 
for  the  central  deflection  was  obtained  by  Ref.  [10]  as 
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The  following  numerical  values  are  taken: 

PQ  = 1000  lb,  t » 20  x 10‘4 * 6  second,  E^  = 30  x 106  psi  and 
the  beam  is  assumed  to  be  j - in  x ^ - in  in  cross-section  and  30  in. 
in  length. 

The  analytical  results  are  obtained  by  retaining  1000  terms  in 
Eqs.  (18)  and  (19).  A check  shows  that  there  is  virtually  no  difference 
between  a 250-term  solution  and  a 1000-term  solution.  Accordingly,  we 
can  practically  take  the  1000-term  solution  as  the  exact  solution. 

Due  to  the  symmetry,  the  higher  order  beam  element  model  is  generated 
for  only  one-half  of  the  beam.  The  solution  is  carried  out  for  N = 5, 

10,  20,  and  50,  where  N is  the  number  of  elements  used  for  the  half  beam. 

For  the  purpose  of  comparison,  the  solution  of  the  conventional  4x4 
element  [11]  is  also  carried  out.  The  central  deflections  at  various 
times  are  shown  in  Table  1.  It  can  be  seen  that  the  higher  order  element 
is  superior  to  the  4 x 4 element.  To  manifest  the  far  more  superiority 
of  the  higher  order  element  in  obtaining  stresses,  the  bending  stress  at 
the  central  span  of  the  beam  is  carried  out  and  shown  In  Table  2.  It  is 
easy  to  recognize  that  the  accuracy  of  the  20-higher-order-element  solution 
is  about  the  same  as  that  of  the  50-  4 x 4 element.  It  takes  CDC  6500 
105  seconds  central  processor  time  for  the  former  and  320  seconds  for 
the  latter.  Thus,  for  the  same  degree  of  accuracy,  the  higher  order  element 
model  uses  about  one-third  the  time  as  would  be  used  by  the  4 x 4 element 
model . 


4 . Elastic  Impact 

When  subject  to  the  Impact  of  a mass,  the  beam  receives  an  impulsive 

force  which  is  the  contact  force  between  the  mass  and  the  beam.  The 
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Table  2.  Bending  Stress  (psi)  at  the  Mid-span  of  the  Beam 
of  an  Isotropic  Material. 


calculation  of  the  contact  force  depends  crucially  on  the  local  deformation 
at  the  contact  region.  The  local  deformation  in  turn  is  affected  by  the 
deflection  of  the  beam.  Timoshenko  £12]  Initiated  a basic  approach  in  this 
area  by  combining  the  Hertzian  contact  force  with  the  Bernoulli -Euler 
beam  theory  to  analyze  the  transient  dynamic  response  of  a beam  subjected 
to  the  impact  of  an  elastic  sphere.  Several  investigators  later  continued 
Timoshenko's  work  by  either  simplifying  the  computational  procedure  or 
employing  the  Timoshenko  beam  theory  in  the  analysis  [13-14].  Many  ex- 
perimental works  Jiave  also  been  conducted  [15-16].  In  all  the  previous  works, 
it  was  assumed  that  the  vibration  of  the  striking  mass  could  be  neglected. 
This  assumption  Is  supported  by  the  result  obtained  by  Rayleigh  that  very 
little  vibration  is  induced  in  an  oscillating  system  under  the  influence 
of  forces  of  duration  long  in  comparison  with  its  natural  periods. 

Consider  an  elastic  sphere  striking  the  beam  in  the  transverse 
direction.  If  the  velocity  of  the  sphere  is  not  too  high,  the  local  de- 
formation between  the  sphere  and  the  beam  can  be  regarded  elastic,  and, 
as  a consequence,  the  Hertzian  law  of  contact  can  be  employed.  Let  F be 
the  contact  force  and  a be  the  relative  approach,  then  the  Hertzian  law 
Is  given  by 

F = k2a3/2  (21) 

where  the  contact  point  spring  constant  kg  is  given  by 


4 
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1/2 


l-vh2 
+ f — — ) 


(22) 
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In  Eq.  (22),  R$  is  the  radius  of  the  sphere;  v$,  E$  and  vb,  Eb  are  the 
Poisson's  ratios  and  the  Young  s moduli  of  the  sphere  and  the  beam,  re- 
spectively. A derivation  of  Eq.  (21)  can  be  found  in  Ref.  [17).  The 
relative  approach  a can  be  expressed  in  terms  of  the  displacements  of  the 
mass  and  the  beam  as 


a = w - v(xQ)  (23) 

where  w is  the  displacement  of  the  sphere  measured  from  the  position  of  the 

initial  contact  and  v(xQ)  is  the  displacement  of  the  beam  at  the  point  of 

contact  x = x„. 

o 

The  equation  of  motion  governing  the  mass  is  given  by 

m$w  = -k2[w  - v(xQJ3/2  (24) 

Thus,  the  motion  of  the  sphere  is  nonlinearly  coupled  with  that  of  the 
beam. 

From  Eq.  (24)  we  have  at  time  t + At 

VW  = - Vat(xo,:>3/2 

and 

Ft+At  * " msWt+At 

where  Ft  t denotes  the  contact  force  at  t+At.  As  a first  estimate,  we 
assume  that  the  values  of  *t+At  and  vt+At(xQ)  can  be  approximated  by 

Wt+At  * *t  +At*t  + \ wt  (27) 


(25) 

(26) 
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and 


Vt+At  (V  = Vt  U0J  + (X0}  + \ ^ t (X0}  (28) 


respectively.  Substituting  Eqs.  (27)  and  (28)  in  Eq.  (25),  we  obtain 
H*.  . in  terms  of  the  known  results  at  the  previous  time  t.  The  contact 

t+At 

force  Ft+  t can  subsequently  be  obtained  from  Eq.  (26).  Using  this 

expression  for  the  contact  force  in  the  established  beam  program,  the 

displacement,  velocity  and  acceleration  at  each  nodal  point  can  be  computed. 

However,  more  accurate  values  of  w.A  . and  v4.J.1*(x„)  can  be  obtained 

t+At  t+At  0 

from  the  first  approximate  values.  Using  the  formula  given  by  Eq.  (10), 
we  can  write 

Wt+At  = Wt  + At  “t  + I At2vit  + l At2  "t+At  {29) 

Vt+At  (X0}  = Vt  (X0}  + At  ;t(X0}  + I At2^  (X0}  + i At2''t+At{X0)  (30) 


where  w*.  . and  v. are  obtained  from  the  first  round  solution. 

t+At  t+At 


Sub- 


stituting Eqs.  (29)  and  (30)  in  Eq.  (25)  and  subsequently  in  Eq.  (26)  we 
obtain  a modified  contact  force  Ft+At*  which,  in  turn,  is  used  to  yield 
a more  accurate  solution  for  the  beam.  Such  iterative  procedure  can  be 
carried  on  till  the  convergence  criterion  is  met. 


5.  The  Timoshenko  Problem 

In  1913,  Timoshenko  [12]  investigated  a 0.394  - in  x 0.394  - in  x 
6.05  - In  (1-cm  x 1-cm  x 15.35  cm)  simply-supported  beam  subjected  to  the 
transverse  Impact  of  a steel  ball  of  0.394  In  radius.  The  Initial 
velocity  of  the  sphere  was  0.394  In/sec.  The  material  constants  taken 
were 
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(31) 


p = 0.0089  slug/in3,  Efa  = 31.3  x 106  pst,  v = 0.289 

In  Timoshenko's  solution,  the  resulting  nonlinear  integral  equations  were 

solved  by  a stepwise  numerical  integration  procedure. 

Due  to  the  symmetry,  only  half  of  the  beam  is  taken  in  the  present 

.8 

analysis.  The  time  step  At  = 5 x 10  second  is  used.  The  contact  force 
history,  the  central  deflection  of  the  beam  and  the  displacement  of  the  sphere 
are  computed  till  the  contact  ceases.  The  results  are  shown  in  Fig.  1. 

It  is  found  that  the  present  solutions  aqree  very  well  with  those  obtained 
by  Timoshenko.  The  rebound  velocity  of  the  sphere  at  the  end  of  the 
contact  is  -0.1192  in/sec  as  compared  with  -0.1198  in/sec  obtained  by 
Timoshenko. 

A similar  problem  can  be  found  in  Ref.  [17],  where  a central  transverse 
impact  of  a 1 /2-in-diameter  steel  ball  on  a 1/2-in  x 1/2-in  x 30-in  simply- 
supported  steel  beam  was  studied.  The  initial  velocity  of  the  ball  was 
assumed  to  be  150  ft/sec.  The  material  constants  were 

E$  = E^  = 30  x 106  psi,  = 0.3,  p$  = = 0.0088  slug/in3  (32) 

In  Fig.  2,  the  contact  force  obtained  by  using  65  harmonics  and  time 

increment  At  = O.lu  sec  Is  reproduced  from  [17]. 

The  solution  according  to  the  present  finite  element  Is  obtained  by 

employing  50  elements  for  half  of  the  beam.  The  time  step  At  is  chosen 
.8 

to  be  5 x 10  sec.  The  results  are  also  shown  In  the  figure  for  the 
sake  of  comparison.  The  agreement  between  these  two  solutions  Is  excellent 
over  a wide  range  of  time.  However  the  maximum  force  obtained  from  the 
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CONTACT  FORCE  , LB. 


Fig.  1.  Contact  Force  and  Displacements  of  the  Timoshenko  Problem 
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Fig.  2.  Contact  Force  and  Displacements  of  the  Goldsmith  Problem 
(with  ."t  = 0.1  sec  and  65  Harmonics). 
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finite  element  method  is  about  2A  higher  than  that  obtained  by  [J7].  It 
was  revealed  in  [17]  that  if  At  was  taken  to  be  2.735;i  sec,  then  the 
maximum  contact  force  would  be  reduced  by  about  9r . Thus,  the  foregoing 
discrepancy  could  have  been  reduced  if  smaller  time  steps  were  used  by 
Gc ldsmi th . 

6.  Impact  with  Permanent  Indentation 

Permanent  craters  can  be  generated  even  at  relatively  low  impact 
speeds.  According  to  Davis  [is],  a permanent  crater  formed  when  a 3/8- 
in-diameter  steel  sphere  dropped  from  a height  of  0.15  inch  onto  an  armor 
plate  which  had  a Grinell  hardness  321.  Therefore,  it  is  more  realistic 
to  use  a dynamic  plastic  force-indentation  law  instead  of  the  Hertzian 
law  of  contact  for  the  impact  problem. 

One  of  the  approaches  in  obtaining  the  dynamic  plastic  force- 
indentation  law  was  proposed  by  Barnhart  and  Goldsmith  [19].  A sphere 
was  dropped  onto  a steel  block  from  heights  of  1 to  13  feet.  Initial 
velocity,  rebound  velocity  and  permanent  crater  depth  were  measured  in  each 
test.  It  was  assumed  that  ( L;  the  initial  kinetic  energy  was  completely 
transformed  into  the  kinetic  energy  of  rebound  and  the  energy  needed  to 
form  the  permanent  crater,  ( 2 ) the  rebound  energy  of  the  sphere  was 
derived  solely  from  the  elastic  strain  energy  stored  equally  in  the  sphere 
and  the  block  during  the  indentation  process,  ( 3 ) the  force  law  could  be 
derived  from  the  modified  Hertzian  law  of  contact:  for  indentation  process, 


F = Nu 


and  for  recovery  process. 


u - u 0 

F ■ HH 

m r 


(33) 
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where  subscript  m denoted  the  state  of  maximum  relative  aoproach,  a 

was  the  permanent  crater  depth,  q could  be  taken  as  3/2  as  in  the  Hertzian 

law,  and  N and  n were  determined  from  iterations  of  the  following  equations: 


N n+1 
n+T  ‘ m 


1 ms  “o  ■ 1 ms  n 


(35) 


1 „ *2  _ m / 

j n,c  («„ 


In  Eqs.  (35)  and  (36),  ms  was  the  mass  of  the  sphere,  wQ  and  w^  were  the 
initial  and  rebound  velocities  of  the  sphere,  respectively. 

Using  the  force-deformation  laws  given  by  Eqs.  (33)  and  (34),  we 
consider  a 1/2-in-diameter  steel  ball  striking  at  the  mid-span  of  a simply- 
supported  beam.  The  steel  beam  has  the  dimensions  1/2-in  x 1/2-in  x 30-in 
and  the  material  constants  are  given  by  Eq.  (32).  The  initial  velocity  of 
the  sphere  is  150  ft/ sec. 

From  Ref.  [ 17]  we  found  that  the  F-u  law  assumes  the  form 


F = 1.29  x 10f-  a1  -128  0 - a - « 

m 


(37) 


and  from  Ref.  [16],  we  found  that  the  average  crater  depth  was 
ur  = 0.0123  in.  Thus,  Eq.  (34)  takes  the  form 


am  - a - 0.0123 
m 


(38) 


The  computer  program  developed  earlier  for  the  elastic  impact  can 
be  readily  extended  to  the  present  case.  During  the  indentation  period. 
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the  Hertzian  law  of  contact  used  in  the  elastic  impact  is  simply  replaced 
by  Eq.  (37),  and  during  the  recovery  process  Eq.  (38)  is  used.  It  should 
be  noted  that  Fm  and  am  appearing  in  Eq.  (38)  can  be  obtained  from  the 
results  for  the  identation  process  when  the  maximum  relative  approach 
is  reached. 

The  results  obtained  based  on  50  finite  elements  and  At  = 0.05y 
second  are  shown  in  Fig.  3.  The  contact  force  given  by  Goldsmith  [17]  is 
also  plotted  for  comparison.  Again,  the  agreement  is  good  over  the  whole 
period  except  in  the  vicinity  of  the  maximum  contact  force.  This  dis- 
crepancy can  be  explained  in  the  same  manner  as  in  the  elastic  case 
described  in  Section  II. 5. 


7.  Finite  Element  for  a Laminated  Composite 

If  the  Bernoull i -Euler  beam  theory  is  employed  to  model  the  symmetric 
cross-ply  laminated  composite,  then  its  stiffness  matrix  can  be  obtained 
from  Equation  (7)  by  replacing  E^I  by  the  bending  stiffness  of  the  laminate 
given  by  h/? 


Db  = 


bE.j  z2  dz 


-h/2 


(39) 


where  b is  the  width  of  beam,  E^  is  the  Young's  modulus  of  the  ith  lamina 
in  the  x-direction,  z is  the  thickness  coordinate,  and  h is  the  thickness 
of  beam.  It  is  noted  that  E^  = EL(the  longitudinal  Young's  modulus)  if  the 
fibers  direction  in  the  ■ th  lamina  is  parallel  to  the  x~axis.  If  the  fi- 
bers are  transverse  to  x-axis,  then  E..  = Ey  ( the  transverse  Young's 
modulus) . 

The  indentation  law  for  an  elastic  isotropic  sphere  on  the  laminated 
composite  is  assumed  to  have  the  same  expression  given  by  Eq.  (21).  Such 
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Fig.  3.  Contact  Force  and  Displacements  of  a Simply-Supported  Beam 
with  Permanent  Indentations  after  Impact. 
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formula  was  found  to  be  valid  by  Willis  [22]  for  a rigid  sphere  pressed  on 
a transversely  isotropic  halfspace.  In  the  present  work,  the  spring 
constant  is  taken  as 


4 

3 


(40) 
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SECTION  III 


IMPACT  RESPONSE  OF  COMPOSITE  MATERIALS 

We  now  proceed  to  use  the  finite  element  and  the  modified  Hertzian 
contact  law  to  study  the  dynamic  response  of  a laminated  beam  and  the 
energy  imparted  from  the  striking  mass  to  the  beam.  It  is  clear  that 
our  main  objective  is  to  predict  the  work  done  by  the  projectile  on  the 
beam.  It  should  be  noted  that  the  total  amount  of  energy  of  the  mass 
is,  in  general,  absorbed  in  two  forms,  namely  the  vibrational  energy 
stored  in  the  form  of  vibration  of  the  laminate  and  the  energy  dissipated 
in  the  contact  zone.  For  a hard  projectile,  damage  generally  occurs 
in  the  contact  area  and  the  vibrational  energy  which,  in  this  case,  does 
little  damage  to  the  beam  will  be  damped  out.  In  view  of  this,  it  is 
desirable  to  separate  these  two  types  of  energies  received  from  the 
striking  mass  and  relate  the  "damage  energy"  to  the  reduction  of  strength 
of  the  laminate. 

1.  The  Elastic  Response. 

Consider  cantilever  beams  of  laminated  composites  subjected  to  the 
impact  of  a mass  at  the  midpoint  of  the  length.  By  using  the  finite 
element  method  described  in  Section  II,  the  entire  response  history  of  the 
beam  and  mass  can  be  obtained.  In  Figs.  4-14,  results  for  various  approach 
velocities  are  presented  for  Scotch  ply/1002  laminated  composites.  The 
elastic  constants  are 

El  = 5.7  x 106  psi,  Ey  = 1.7  x 106  psi,  p = 0.00202  slug/ in } 

The  projectile  is  a steel  ball  of  diameter  0.177  in.  Twenty-four 
elements  are  used  with  At  = 0.1  u second.  The  same  dimensions  of 
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Mg.  4.  Impact  Response  History  of  the  Laminated  Cantilever  of 
0.5"  x 0.097"  x 6"  for  Approach  Velocity  of  150  fps. 
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Fig.  5.  Impact  Response  History  of  the  Laminated  Cantilever  of 
0.5"  x 0.097"  x 6"  for  Approach  Velocity  290  fps. 
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Fig.  9.  Impact  Response  History  of  the  Laminated  Cantilever  of 
1"  x 0.102“  x 6"  for  Approach  Velocity  775  fps. 
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Fig.  13.  impact  Response  History  of  the  Laminated  Cantilever  of 
1.5"  x 0.103"  x 6"  for  Approach  Velocity  880  fps. 
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Fig.  14.  Impact  Response  History  of  the  Laminated  Cantilever  of 
1.5"  x 0.104"  x 6"  for  Approach  Velocity  904  fps. 
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specimens  and  materials  were  also  used  by  Husman  et  al  in  their  experiments 
[5]. 


The  energy  imparted  from  the  sphere  to  the  laminate  during  the  period 
of  contact  can  be  evaluated  in  several  ways.  We  can  either  compare  the 
initial  and  rebound  velocities  of  the  mass  or  the  contact  force  and  dis- 
placement histories.  We  have 


or 


KE  = \ ms(w02  ~ «f) 


K£  = 


F dw 


(42) 

(43) 


where  w^  is  the  displacement  of  the  mass  when  the  contact  ceases. 

Of  course,  one  cannot  expect  these  elastic  solutions  to  compare 
favorably  with  the  experimental  results.  The  reason  is  obvious:  The 
inelastic  response  in  the  contact  zone  is  too  substantial  to  neglect. 


2.  Estimate  of  the  Imparted  Energy 

The  finite  element  procedure  can  be  used  for  any  kind  of  inelastic 
contact  law.  In  the  future,  if  such  laws  are  established,  the  finite 
element  solution  should  be  able  to  yield  more  accurate  response  histories. 
However,  for  the  FOD  of  composites  problems,  we  are  mainly  concerned 
with  the  energy  absorbed  by  the  laminate  rather  than  Its  dynamic  response 
history.  With  this  In  mind,  we  are  able  to  utilize  the  special  material 
properties  of  the  fiber-reinforced  composite  to  estimate  this  energy  with 
surprisingly  accurate  results. 

From  all  the  available  experimental  data,  the  stress-strain  relations 
of  fiber-composite  composites  behave  almost  linearly  up  to  failure.  Only 


35 


the  longitudinal  shear  exhibits  somewhat  nonlinear  behavior.  Since  in  the 
present  consideration  the  longitudinal  shear  deformation  is  negligible,  the 
linear  relations  are  valid  in  the  loading  stage.  Only  when  the  indentation 
reaches  its  maximum  depth  and  starts  to  recover,  the  linearly  elastic 
relations  should  be  subject  to  reexamination.  If  we  assume  that  the  elastic 
recovery  of  deformation  is  negligible  as  in  rigid  plasticity,  then  the 
relation  of  the  contact  force  and  the  indentation  a can  be  represented  by 
the  curve  shown  in  Fig.  15.  In  other  words,  the  contact  force  will  drop 
to  zero  after  it  passes  the  maximum  value  of  a in  the  elastic  solution;  and 
in  the  subsequent  time,  the  mass  does  no  work  to  the  beam.  The  total 
amount  of  work  done  by  the  projectile  becomes 


KE 


T 


f'W  Fdw 


(44) 


where  wmax  is  the  displacement  of  the  sphere  at  F = F . It  should  be 
noted  that  wmax  is  not  the  maximum  value  of  w. 

Using  Eq.  (44),  we  compute  KEy  and  tabulate  the  results  in  Table  3. 

Our  analytical  results  are  excellent  as  compared  with  the  experimental 
results  obtained  by  Husman  et  al  [5]. 

3.  The  Damage  Energy 

Since  not  all  the  work  done  by  the  projectile  can  be  converted  to 
energy  dissipated  In  the  damage  zone,  the  total  Imparted  energy  KEy 
cannot  be  used  to  account  for  the  amount  of  damage  received  by  the  laminate. 
A more  pertinent  quantity  is  the  damage  energy  defined  by 


KE 


D 


F da 


(45) 
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al  Imparted  and  Damage  Energies 
Laminated  Beams  of  a Glass-Epoxy 
Composite. 


The  rest  of  the  work  done  i&  stored  in  the  form  of  vibrational  energy 
given  by 

KEy  = KET  - KED  = f V,nax  F dv  (46) 

' o 

where  v is  the  displacement  of  the  beam  at  the  contact  point  when 

F = F . The  results  are  also  shown  in  Table  3. 

max 


I 


SECTION  IV 


DISCUSSIONS  AND  RECOMMENDATIONS 

Damages  of  turbine  fan  blades  of  advanced  composites  due  to  impact 
of  foreign  objects  are  to  be  related  to  their  residual  strengths.  In  this 
report,  an  effort  is  made  to  develop  a finite  element  method  for  determin- 
ing the  amount  of  energy  that  causes  damages  in  the  area  of  impact.  The 
model  is  derived  based  upon  linear  stress-strain  relations  and  a modified 
Hertzian  contact  law.  The  results  are  excellent  in  comparison  with 
existing  experimental  data  for  glass-epoxy  laminated  composites.  With 
this  analytical  procedure,  more  general  type  of  laminates  can  be  con- 
sidered and  optimal  design  with  respect  to  impact  strength  of  blades  is 
feasible.  The  greatest  advantage  of  this  model  resides  in  its  simplicity 
and  without  using  inelastic  material  behaviors  which  have  been  deemed 
necessary  in  solving  FOD  problems. 

However,  in  adopting  this  model,  several  restrictions  must  be 
imposed.  First,  the  blade  can  only  be  modeled  as  a beam.  Second,  the 
projectile  should  be  relatively  hard  so  that  the  damage  zone  is  confined 
in  the  vicinity  of  the  contact  region.  The  applications  are  more  suitable 
for  impact  of  hailstones  and  gravels  rather  than  birds.  The  response  of 
blades  to  impact  of  birds  is  substantially  different  from  that  of  the 
hard  objects.  A soft  body  such  as  a bird  usually  leads  to  a longer 
contact  time  and  is  able  to  produce  larger  amplitude  stress  waves  in  the 
blade.  Such  waves  with  high  stress  levels  could  produce  failure  of  the 
blade  at  distant  regions  such  as  the  root,  and  the  local  damage  might  be 
less  significant.  Furthermore,  the  bird  itself  could  absorb  a substantial 
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amount  of  energy  during  its  destruction.  In  using  fi.iite  element  for  soft 
body  impact,  the  finite  element  for  the  blade  must  be  modeled  to  account 
for  large  deflection  and  the  soft  body  be  modeled  by  viscoplastic  finite 
elements.  This  part  of  the  work  is  underway  and  will  be  reported  in  the 
future. 

In  order  to  model  more  realistically  an  actual  turbine  fan  blade,  a 
plate  finite  element  should  be  used.  The  two-dimensional  plate  element 
can  describe  the  twist  geometry  of  the  fan  blade  by  using  the  usual  co- 
ordinate transformation  for  finite  elements. 

Finally,  although  it  is  of  less  practical  interest,  the  elastic 
unloading  can  be  incorporated  in  the  present  analysis  in  a crude  manner. 

As  mentioned  previously,  unloading  occurs  after  the  contact  force  passes 
its  maximum  value  as  shown  schematically  in  Fig.  16.  Exact  unloading 
path  is  unclear  at  this  stage.  We  propose  to  estimate  the  elastic  strain 
to  be  recovered  as 


ee  = Y/Et  (47) 

where  Y is  the  transverse  strength  of  the  fiber-reinforced  composite. 

The  total  recovery  displacement  at  the  is  given  by 

°e  = \ h Ee  (48) 

The  one-half  factor  is  used  to  account  for  the  fact  that  the  compressive 

strain  over  the  thickness  of  beam  Is  not  uniform  and  the  average  value 

is  taken.  With  this  consideration,  the  projectile  remains  in  contact  with 

the  beam  till  the  value  of  a passes  a and  reaches  <*„  -ct  . Thus,  the 

r max  max  e 
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total  work  done  by  the  projectile  becomes 


KE* 


max 


F dw  + 


1 p w 

2 max  we 


(49) 


where  wg  is  the  displacement  of  the  projectile  corresponding  to 

a = “max  " “e  after  Passin9  « = “max-  It  is  noted  that  one-half  is 
added  to  the  second  term  to  approximately  account  for  the  fact  that  the 
contact  force  should  drop  from  Fmax  to  zero  in  the  region  of  unloading. 
The  damage  energy  KE^  can  be  corrected  in  a similar  manner.  We  have 


KV  = KE0  - 7 Fmax  “e 


(50) 


In  the  calculation,  we  choose  Y = 4 x 103  psi.  Both  KET*  and  KEQ*  are 
presented  in  Table  3.  Some  improvement  is  noted. 
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